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The abundance of a species’ population in an ecosystem is rarely stationary, 
often exhibiting large fluctuations over time. Using historical data on marine 
species, we show that the year-to-year fluctuations of population growth rate 
obey a well-defined double-exponential (Laplace) distribution. This striking 
regularity allows us to devise a stochastic model despite seemingly irregular 
variations in population abundances. The model identifies the effect of re¬ 
duced growth at low population density as a key factor missed in current ap¬ 
proaches of population variability analysis and without which extinction risks 
are severely underestimated. The model also allows us to separate the effect of 
demographic stochasticity and show that single-species growth rates are dom¬ 
inantly determined by stochasticity common to all species. This dominance— 
and the implications it has for interspecies correlations, including co-extinctions 
—emphasizes the need of ecosystem-level management approaches to reduce 
the extinction risk of the individual species themselves. 

Keywords: complex systems, time series, population dynamics, growth rate statistics, stochas¬ 
tic processes, nonlinear dynamics 
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1 Introduction 


Assessment of extinction risk and biodiversity loss is a central problem in ecology, which has 
direct implications for ecosystem management practices 0 and policies for the exploitation of 
natural resources [2j. It is estimated that currently up to 0.1% of all known species go extinct 
every year, which is over one thousand times above the background extinction rate observed 
in fossil records [3j. Whether caused by habitat degradation, interspecies competition, climate 
change, overexploitation, or the introduction of exotic species, the majority of such extinction 
events have not been anticipated. Existing inventories of the global conservation status, such 
as the IUCN Red List []4]|, are believed to include only a fraction of all endangered species 
and, even for those, significant uncertainty remains on their actual extinction risk. The central 
difficulty is that wild populations generally do not exist in a steady state from which minute 
deviations or trends could be detected. Instead, they tend to exhibit fluctuations over time (5]-{9j. 

Temporal variations in population abundance may be caused by species interactions, en¬ 
vironmental changes, migration patterns, intrinsic nonlinearities, and/or human exploitation 
pO[jT2| , and are sometimes sufficiently irregular to be regarded as stochastic. Such irregu¬ 
lar time dependence, when combined with unavoidably imperfect sampling of the population, 
poses considerable challenges for population forecast. Significant previous research has focused 
on determining the frequency composition (so called “noise colour”) of such fluctuations and 
their correlations with the environment [13] 141. Despite their immediate implications for sus¬ 
tainable ecosystem exploitation and management, much less understanding has been generated 
about the factors that influence the growth rate of individual species. 


Perhaps not surprisingly, historically there has been a large number of both false positives 
and false negatives in the assignment of extinction status and extinction risk. For instance, the 
New Zealand’s bird takahe, which was considered extinct by the end of the 19th century, was 
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rediscovered in the wild 50 years later ©• This is one of now many known examples of a 


Lazarus taxon [16], in which a sparse population passed undetected for an extended period of 
time. The passenger pigeon in North America, on the other hand, went from an abundant pop¬ 
ulation to functional extinction in less than 20 years—which then led to actual extinction in the 


early 20th century [17]. Taken together, the picture that emerges is one in which the survival of 


a species appears to depend very subtly on the species’ own abundance [18]. This picture is fur¬ 


ther complicated by the possibility of co-extinctions [19j20|, whose actual role beyond directly 


dependent species (such as predator-prey and parasite-host) remains elusive [21-j23|. Compared 
to the terrestrial case, extinctions of marine species caused by habitat loss or human exploita¬ 
tion have been rare, in part due to stabilizing effects such as changes in fish catchability at low 
population densities. Yet, the pace of marine defaunation is likely to accelerate dramatically as 


the strength and scope of human impact on marine ecosystems grow [24]. This underscores the 
need for a quantitative understanding of the growth dynamics of marine populations, including 
extinction risks. 


2 Results 


Here, we examine the dependence of growth rate on population abundance and stochastic fac¬ 
tors, both when only the temporal variability of individual species is observed and when the 
dynamics of the entire ecosystem is taken into consideration. We base our analysis on data of the 
Large Marine Ecosystems (LME) portion of the Sea Around Us Project (http://www.seaaroundus.org) 


[25], which is a global-scale database on species abundance in 66 marine ecosystems. For each 
ecosystem, the data consist of the annual quantities (in tonnes) of its 12 most abundant species 
caught by fisheries over the period 1950 to 2006. Such landing data is a widely used tool for 


making inferences about marine populations [26-29], even if such use is occasionally con¬ 


troversial [30], as factors other than population abundance can influence commercial catches. 
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We will show, however, that our results hold true for marine stock assessments, which inte¬ 
grate research surveys and catch data with independent information (such as mortality rates and 
population size/age structure) to obtain a more accurate estimate of population abundance (sup¬ 
plementary material, Analysis of Stock Assessments). Therefore, we present our analysis based 
on the more widely-available landing data |25|, followed by supplementary validation using 
stock assessment data. 

For each species * in a given ecosystem, we assume that the annual abundance in year t 
can be approximated—up to a scaling factor—by the reported landing, denoted by x t , al¬ 
though we will show that our results do not depend critically on this assumption (supplemen¬ 
tary material, section 1, figures SI-S3). Our object of study is the year-to-year growth rate, 
Tf'- 1 = In [x^/x^. To avoid ill-defined log functions associated with zero landings, we add 
1 to each data point, so that the minimum value of x f 1 is 1. Without any further assumptions, 
the year-to-year change in the species’ population abundance can always be written as 


x 


(0 
t +i 



( 1 ) 


where = r (l) + represents the growth rate in year f decomposed into the average 

fW and standard deviation crW 0 f the growth rate (calculated over the entire time series) and 
a time-dependent factor Q at time t. The term Q is thus the normalised growth rate fluctu¬ 
ation. For example, Eq. [l] reduces to the classical Ricker model [31j if is a deterministic 
decreasing linear function of the abundance which, as shown below, does not hold true for 
the marine ecosystems we consider. In contrast with previous studies, here we will make no a 
priori assumptions on the growth rate, instead deriving its properties directly from the data. 

Figure [ljpresents empirical properties of x/f’ for all 66 ecosystems we consider. We first note 
that while the autocorrelation of the (log) population abundances In xf* is significant and close 
to one for successive years, as expected (figure [l]4), the autocorrelation of the corresponding 
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growth rates is comparatively low, with a majority of species having autocorrelation less than 0.2 
in magnitude. This surprising observation indicates that the time-dependent component of the 
growth rate, can be regarded as a random variable drawn from an appropriate distribution 
that is nearly stationary. The data show that, to an excellent approximation, this distribution is 
given by a double-exponential function 

F i(0 = ^exp(-v / 2^|), (2) 

also known as the Laplace distribution (figure [T^-C). This is itself an important, novel finding, 
which further allows us to devise a statistical model, as discussed below. We have verified that 
our addition of l’s to the zero-landing years does not affect this distribution, which incidentally 
also governs the fluctuations of growth rate for the population abundance of each ecosystem as 
a whole (figure [lp). In addition, we have rigorously confirmed the fit of the individual species’ 
growth fluctuations to the Laplace distribution versus the null hypothesis of a normal distri¬ 
bution, according to standard goodness-of-fit tests such as the Kolmogorov-Smirnov test and 
Akaike information criterion. These tests show that the former distribution is a significantly 
more plausible explanation than the latter distribution for the growth rates of a large majority of 
species (supplementary material, figure S4). 

Having established that the normalised growth rate fluctuations obey a stationary stochastic 
process described by a parameter-free Laplace distribution, we can propose Eqs.[l]j2]themselves 
as a stochastic model of population abundance. This model is expected to be appropriate away 
from the floor abundance x) = 1. However, the growth rate at the floor is much smaller than 
at any other abundance, with probability 89% of being zero, as shown in the inset of figure [Ip. 
Naturally, because abundance is measured in integer units of landing tonnes, x t = 1 only 
indicates that the population is very low but not necessarily that it is zero. For this reason (and 
possibly due to migration and sample biases), the apparent “extinctions” considered here are 
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local in space and in time. Indeed, the populations generally recover to detectable abundances 
in the systems under consideration. However, they do so more slowly than predicted by the 
overall growth rates. This reduced growth rate at low population density is consistent with the 
Allee effect (or depensatory population dynamics), which is a scenario previously observed for a 
number of species in diverse ecosystems fl32| . Further analysis would be needed to conclusively 
address that particular issue in this case, which falls outside the scope of this paper. It should 
be noted, nevertheless, that the apparent reduced growth rate at low population abundances is 


counter to most existing models, including the classical Ricker model [33]. A recent study of 
exploited marine species stocks included in the RAM Legacy Stock Assessment Database 0 
indicates that this missing element is in fact the probable explanation for the slow recovery 
of depleted stocks when compared with predictions from models commonly used in fisheries 


management [35J. 

We incorporate the floor effect into our model by taking the probability distribution of 
to be 

>i(0, 


(0 


HO = 


for Xf 1 ' 1 > 1, 


(1 — p^) x 5(£ — _|_ pW x 2Fi(£), fora;)*' 1 = 1 and £ > 0, 


(3) 


and _F 2 (£) = 0 otherwise. The parameter pW i s a measure of the recovery probability once the 
species reaches the floor abundance, 6 is the Dirac delta function, and ^ = —fW/cr^ is de¬ 
fined so that = 0 when ^ l >) = H (for simplicity in equation (3), we used the approximation 
> 0 instead of the exact condition since fW—and hence ^—is typically close 

to zero). Here, the average growth rate fW and standard deviation crW are estimated for x[ l) > 1 
and the recovery probability pW i s estimated for x, !> = 1 ( Methods , Parameter estimation). 
The model defined by Eqs.[[|{3] offers projections on future population abundance based on past 
abundance, which in turn can be used to assess risk of pseudo-extinctions , defined as crossings 
below a given fraction 6 of the historical maximum population. 
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Figure [2] validates our model against empirical data {Methods, Model predictions and vali¬ 
dation). It shows, in particular, that the floor effect is crucial for the excellent agreement found 
between the pseudo-extinction predicted and the ones actually observed over the same period 
(figure [2j\). We have tested that the heightened agreement with the empirical data holds re¬ 
gardless of the exact value used for the floor abundance (which also represents the limit of 
resolution in the data), insofar as it is not too large (figure[3]). This is significant because current 
approaches for population variability analysis, including specialised commercial software, gen¬ 
erally do not account for the exceptional case of growth at very low population densities (even 


though minimal viable population sizes are often assumed [36]). Neglecting the floor effect not 
only mis-predicts pseudo-extinctions observed in the past, but also tends to severely underes¬ 
timate the risks of future pseudo-extinctions (figure [2|3). For example, the popular deepwater 
redfish {Sebastes mentella) in Baffin Bay appears to have recovered from very low population 
abundances, but historical data indicate that when the abundance of this species approaches 
zero it has a high probability of remaining extremely low for extended periods (figure [2p). 
Neglecting this reduced growth leads to a significant underestimation of the pseudo-extinction 
risk. This underestimation is even more pronounced for a number of other species, such as the 
flatfishes (Pleuronectiformes spp.) in the Antarctic (figure [2p)). Interestingly, for some species 
neglecting the floor effect may actually lead to an overestimation of the pseudo-extinction risk. 
For example, the blue mussel {Mytilus galloprovincicilis) in the Iberian Coast reached the floor 
abundance in the year 2004, but started recovering immediately (figure|2j3). Because the species 
exhibited positive growth the only time it reached the floor, predictions that take this into ac¬ 
count naturally lead to an estimate of pseudo-extinction risk that is smaller and more reliable 
than predictions that do not. 

Central to our analysis is the fluctuation of the growth rate, modeled as a stochastic term Q . 
But what is the relation between the stochasticity of different species i in the same ecosystem? 
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This question can be addressed by making the ansatz that the Q are not independent but are 
instead drawn from an appropriate joint distribution in which the (marginal) distributions for 
the individual species follow Eq. [3] while having correlation a 2 with those of the other species 
(Methods, Estimation of common stochasticity). We then consider a range of values of a, 
which represents the portion of stochasticity common to all species, and analyse the integrated 
impact on the population abundance of each ecosystem as a whole. To give all the species 
comparable weight, we focus on the average over the individual species’ population weighted by 
the inverse of their average abundance (as in figure [T]D). The occurrence of pseudo-extinctions is 
systematically underestimated when stochasticity is dominantly species-specific, as illustrated 
in figure |4]4 for a = 0. In fact, the agreement of the model with empirical data is significantly 
better when stochasticity is taken to be dominantly common to all species (figure |4j3), and the 
agreement becomes excellent for a = 0.8 (figure [4pk). 

This corroborates the conclusion that variations in growth rate are largely synchronised 
within each ecosystem. We speculate that this synchronisation is partially rooted in external, 
environmental fluctuations, which are known to play a crucial role in the dynamics of terres¬ 
trial populations pTj[T7| . In marine systems, likely environmental fluctuations include the El 
Nino Southern Oscillation, the North Atlantic Oscillation, and riverine flood pulses. For fishes, 
the impacts of these fluctuations can be both direct and indirect—such as via externally-driven 


changes in the lower food web or via ecosystem regime shifts [38]. Other potential sources 
of synchronized growth fluctuations include interspecies interactions as well as correlations 
with human activity—for example shifts in overall fishing effort due to weather or changes in 
regulations. But regardless of the source of the common stochasticity observed here, the di¬ 
rect implication is an increased risk of the otherwise unlikely concurrent collapse of multiple 
species. 
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3 Discussion 


Our findings should be compared with the case studies of the American breeding bird popula¬ 


tions [39] and Hinkley Point’s fish community [40], for which growth rates have been analyzed. 
For both systems the aggregated non-normalised growth rates were found to follow power-law 
distributions, but independent analysis of some of these data has shown that after rescaling (by 
the species standard deviation) to a variable equivalent to the normalised growth rate fluctuation 


the distribution becomes normal [41], which corroborates the conclusion that the growth 


rate distributions of individual species are short-tailed. This is consistent with the previous anal¬ 


ysis of 544 long-term time-series from the global population dynamics database [42], includ¬ 
ing both aquatic and terrestrial populations, which demonstrated that the abundances of most 
species are either lognormally distributed or shorter-tailed than lognormally. If one neglects 
the (significant) year-to-year correlations in abundance, lognormal distributions for individual 
population abundances imply normal distributions for individual species’ growth rates. 

The results presented here, on the other hand, show that normalised growth rate fluctuations 
follow Laplace distributions and this is confirmed to remain true for individual species in the 
ecosystems we consider (supplementary material, figure S4). Had we not normalised the growth 
rates to eliminate heterogeneity across species, the resulting aggregated non-normalised growth 
rate fluctuations would be more fat-tailed than what we observe, without necessarily revealing 
a simple scaling behaviour (supplementary material, figure S5). Importantly, we have verified 
that the observed Laplace distributions for the normalised growth rates are not artefacts of the 
landing data used here as a proxy for population abundance, remaining valid for stock assess¬ 
ment data. This is demonstrated in figure S6 (supplementary material) for the RAM Legacy 
Stock Assessment Database [3J, which is the most complete marine stock assessment catalogue 
available. The estimates of population abundance therein were acquired under controlled set- 
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tings, using a variety of methodologies designed to avoid systematic biases, and in a variety 
of ecosystems, which substantiates the conclusion that Laplace statistics underlie the growth of 
marine populations in general. Note that different types of growth rate distributions are con¬ 
sistent with the lognormally-distributed abundances observed in previous studies; by Eq.|T| the 
(log) population abundance is (up to a constant) the sum of the growth rates in all previous 
years. Thus, by the Central Limit Theorem, any growth rate distribution (including Laplace) 
will eventually lead to normal distributions for the log abundances, provided the growth rates 
are independent and identically distributed with finite variance. However, direct analysis of 
the marine population abundances in this study reveals that their distributions are no closer to 
lognormal than to power laws (supplementary material, figure S7). Altogether, our results are 
significantly different from those suggested by previous studies, and they do not follow from 
existing ecological models. 

Our demonstration that the growth rates of marine species are governed by the Laplace 
distribution, which has a heavier tail than a normal distribution with the same standard devi¬ 
ation, has important implications for the analysis of empirical data. Because the likelihood 
of pronounced fluctuations is larger than expected from normal distributions, large short-term 
fluctuations (over the period of few years) are not necessarily a sign of abnormality as they 
may well be a natural property of the system. This, combined with the observed stickiness to 
the floor abundance, poses additional challenges to the identification of abnormal population 
dynamics. In particular, we have shown that neglecting the floor effect alone already leads to 
substantial underestimation of pseudo-extinction risks. Finally, since Laplace distributions have 


been previously identified in the growth of companies [43], our study establishes a new parallel 


between ecological and socio-economical networks, both of which are characterized by growth 
and competition in the presence of limited resources. As such, our results may also provide 
new insights into common stability mechanisms that govern otherwise disparate systems—as 
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recently proposed in ref. 44 


4 Methods 


Parameter estimation. To estimate the model parameters from a given time series of popula¬ 
tion abundance, we disregard the initial consecutive l’s, if any, as they often correspond to years 
for which no data are available. Using the resulting time series {x t }, we compute the associated 
growth rate time series {ry} from the definition r t = In xt+i — In ay, where we now omit the 
superscript (i) for simplicity (a convention also adopted in the figures). The parameter estima¬ 
tion depends on whether or not we consider the floor effect in the model. For the model without 
floor effect (Eqs.[l]{2]), f and a are simply the sample mean and standard deviation of {rt}. For 

’ we have f =\k\ £teh r t> ° = Y/^piEte/^U-r) 2 , 
and p = 1 — where /, = {t\x t > 1}, I 0 = {t\x t = 1}, / 00 = {t\x t = 1 and x t+1 = 1}, and 

I • I denotes the number of elements in the set. 


the model with floor effect (Eqs. 
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Model predictions and validation. For each time series {ay} of length T, we use its first T 0 
data points as the training set for parameter estimation and the remaining T — T 0 data points for 
validation. We calibrate our model both without and with floor effect, and simulate both variants 
for 1, 000 independent runs. Each run starts with the same initial condition xt 0 and is computed 
for a total of T — T 0 time steps, with ay reset to 1 whenever its estimated value is below 1. For 
both the emprical and simulated data, we calculate the pseudo-extinction risk as the probability 
that x t < OX, where 6 e [^, l) is the pseudo-extinction threshold and X = maxi< 4 <T 0 ay. 
The special case of pseudo-extinction risk at the floor level is defined as the fraction of years 
for which the population abundance is at the floor (i.e., ay = 1). In this study, T = 57 years and 
we choose T 0 = 45 years. 
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Estimation of common stochasticity. We investigate the existence of common stochasticity 
among all species considered in each ecosystem by extending our model as follows. In a given 
year t, for those species that do not remain at the floor abundance (according to Eq. [3]), we 


draw a vector £ t = () of growth fluctuations from a multivariate Laplace distribution [45] 


characterised by (vector) mean 0, (vector) variance 1, and normalised covariance matrix \ . We 
set the diagonal elements of T identically to 1 and the off-diagonal elements identically to a 2 , 
where a is a parameter ranging from 0 to 1. The resulting growth fluctuations for the individual 
species are distributed according to Eq. [i] but the correlation coefficient is a 2 between ^ and 
^ ( t J> for any pair i ^ j. As such, the parameter a can be interpreted as the proportion of common 
stochasticity. The value of a that best represents the empirical data is determined in figure [4] 
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Fig. 1: Statistical properties of population growth rate derived from the Large Marine 
Ecosystems dataset. (A) Distribution of the absolute values of the lag-1 autocorrelation of 
the log-transformed annual population abundance {lnx t } (upper part) and of the growth rate 
{r t } (lower part). For most species, {lna^} exhibits significant autocorrelation whereas {r,} 
does not (B) Normalised growth rate fluctuation for individual species, showing the double¬ 
exponential “tent” shape of the distribution.(C) Same as in panel B, where the distribution 
of positive fluctuations P(£ t |£ t > 0) (solid symbols) and the distribution of the magnitudes 
of negative fluctuations P(—£t\£t < 0) (open symbols) are displayed separately to show the 
symmetry. Both distributions follow an exponential function P(£) ~ exp(—A|£|), with A ~ a/ 2- 
The inset shows the probabilities P(r t = 0| In x t = 0) and P(r t > 0| In x t = 0). When the 
abundance of a species reaches the floor level In x t = 0, it has a high probability (89%) of 
staying at that level the following year, which reveals a strong relation between low population 
and nearly zero growth rate. We note that, while our addition of l’s to the reported populations 
will tend to deflate the apparent magnitude of the growth rate at small population abundances 
when calculated using Eq. |T] the stickiness to the floor reported here is a property of the raw 
data and is not influenced by this transformation. (D) Same as in panel C at the ecosystem level, 
where the ‘population abundance’ of an ecosystem is taken to be the weighted average over all 
of its 12 recorded species, with the weights defined by the inverse of the average abundance of 
the species. The exponent y/2 is within the 95% confidence interval of the linear regressions. 
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Fig. 2: Model validation and predictions for individual-species pseudo-extinction risk. 

(A) Model discrepancy as a function of the pseudo-extinction threshold, 6, for a threshold 
range below 1% of the maximum population. The curves show the average log-distance be¬ 
tween the pseudo-extinction risk in the empirical data and that predicted by the model without 
floor (squares) and the model with floor (stars). The model is parameterized over the period 
1950-1994, and the predictions are shown for the period 1995-2006. The curves correspond 
to averages over all recorded species of the ecosystems under consideration. (B) The risk of 
pseudo-extinction at the floor level, projected for the period 2007-2021 for the model without 
and with floor parameterized over the full period 1950-2006 of available data. The higher inci¬ 
dence of points in the upper right corner above the shaded reason along the diagonal indicates 
that the pseudo-extinction risk is grossly underestimated if the floor is not explicitly accounted 
for (note the logarithmic scale). (C to E) Examples of empirical time series for individual 
species (marked dots in panel B) with zero-landing years leading to large systematic discrep¬ 
ancies between the pseudo-extinction risk predicted by the model without floor (first number) 
and with floor (second number). Pseudo-extinctions are defined as crossings of population 
abundance below the given threshold 6 (measured relative to the historical maximum), while 
pseudo-extinctions at the floor level are operationally defined as the extreme case in which pop- 
ulations reach the floor abundance x) = 1. The associated risk is quantified as the fraction of 
years below threshold or at the floor, respectively. See Methods for details on model validation 
and predictions. 
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Fig. 3: Model validation for alternate floor values and levels of discretization. Counterpart 
to figure [2}\ for different values of the “floor” abundance (measured in tonnes), where dashed 
(continuous) lines indicate the model discrepancy with (without) floor. Incorporating the ob¬ 
served tendency for already low population abundances to remain low significantly increases 
the agreement with the empirically-observed pseudo-extinction risk for low threshold 9, regard¬ 
less of the precise value used to implement the floor effect. For each simulated value of the 
floor, it is assumed that the population abundances are not resolved below the floor level and 
that populations are only resolved to the nearest multiple of this level. 
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Fig. 4: Model validation for ecosystem pseudo-extinction risk. (A) Pseudo-extinction risk 
of the weighted average population of the ecosystem as a function of the threshold 9 : empirical 
data (circles), model for a = 0 (squares), and model for a = 0.8 (stars), where a is the 
portion of the stochasticity common to all species in the corresponding ecosystem. Each curve 
is an average over all ecosystems for the same 45-year parameter-determination period and 12- 
year prediction period used in figure [2]A.. (B) Model discrepancy as a function of or. average 
log-distance (upward triangles) and average signed log-distance (downward triangles) from the 
model curve to the empirical curve for the range of 9 shown in panel A. For comparison we 
highlight the limiting case of a = 1 (dashed line), which is equivalent to applying our model 
directly to the weighted average population of each ecosystem. 
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Regularity Underlies Erratic Population Abundances in Marine Ecosystems 
J. Sun, S.P. Cornelius, J. Janssen, K.A. Gray & A.E. Motter 

Supplementary Material 

1 Effects of Non-Constant Catchability 

For a given species, the potential relationship between the reported catch (landing) in a given 
year, x t , and the underlying population abundance, n t , can be expressed in the general form 

— = (S4) 

n t 

where / is a function that then describes the catchability , or catch per unit population. For 
example, in the classical Type I fishery model, f(n t ) = qE, where q is catchability per unit 
effort and E is total effort [1]. Here we focus on more general functions / and study their ef¬ 
fects on growth rates and population dynamics. We will make the assumption that catchability 
approaches a constant as population size grows large, i.e., f(n t ) —* a as n t —>■ oo. Under this 
assumption, we explore the following three mostly commonly-considered scenarios. 

Baseline. This is the simplest scenario, under which the amount caught is directly proportional 
to the actual population abundance (i.e., f{n t ) = a). Results reported in the main text were 
obtained under this assumption. In this case, the growth rates calculated from catch data are 
equal to the actual population growth rates since x t+ i/x t = n t+ i/n t , and hence the particular 
value of a does not affect the growth rate distribution nor the predictions of pseudo-extinction 
risk made by our stochastic model. 
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Hyperdepletion. In this scenario, catchability decreases at low population abundance, which 
is typically assumed to mean that /(0) =0 and df /dn t > 0 for all n t > 0 [2]. As a concrete 
example, we consider the following functional form: 


/M = 


an t 

n t + b : 


(S5) 


where b > 0 is a parameter. Note that we can use Eq. S4 to solve for population in terms of 
catch in this case as 


n t = 


x t + \Jx\ + 4 abx t 


2 a 


(S6) 


Hyperstability. In this scenario, catchability increases at low population abundance, which is 
typically taken to mean that /(0) > 0 and df /dn t < 0 for all n t > 0 [2]. As a concrete example, 
we consider the following functional form: 


/M 


a + y/a 2 + 4 c/n t 
2 


(S7) 


where c > 0 is an additional parameter. From this we obtain that 


rh = 


x: 


ax t + c 


(S8) 


Figure [ST] illustrates the relationship between catch and underlying population abundance 
for the three scenarios above. In our analysis, we use the parameter values a = 0.1, b = 10 4 , 
and c = 10 2 . These values were chosen based on the typical scales of the catch data in the 
FME dataset. Note that for the purposes of this section, we will calculate growth rates and 









24 


train/validate our model based on the population abundances , which are obtained by transform¬ 
ing the catch data in the LME dataset according to one of the three catchability relations defined 
above. 

Figure [S2] shows the distribution of growth rate fluctuations under the hyperstability and 
hyperdepletion scenarios. In both scenarios, we see that the distribution is indistinguishable 
from the distribution of growth rate fluctuations in the baseline scenario, which we have shown 
to be Laplace (Fig. [I] and Fig. S4 2-D). As such, Laplace statistics should still form the basis 
of our predictive stochastic model under other conceivable catchability scenarios. But how do 
the model’s predictions of pseudo-extinction risk fare when based on population rather than 
reported catch? 

Figure [S3] shows the model-predicted pseudo-extinction risk as a function of pseudo-extinction 
threshold compared to the empirically-observed risk for the hyperstability and hyperdepletion 
scenarios. As for the baseline scenario of constant catchability (Fig. [2]A), the inclusion of 
the floor effect significantly improves our model’s prediction of pseudo-extinction risk at low 
pseudo-extinction thresholds. However, in the case of hyperstability, both models (with and 
without floor effect) systematically underestimate pseudo-extinction risk. 


2 Analysis of Stock Assessments 

We use catch (landing) data for the analysis in the main text because it lends itself to high- 
quality statistics, being available for a large number of years and for a large number of species. 
Indeed, catch data is the only available indicator of the population abundance of most marine 
species. Nonetheless, one must be cautious in using reported landings as a direct proxy for 
population abundance, since there are factors that can affect catches other than changes in un¬ 
derlying population abundance, such as extreme weather events and changes in market demand 
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or fishing effort. 

To address the possibility that these artefacts may have affected our results, we have repeated 
our statistical analysis on the RAM Legacy Stock Assessment Database [3j, which is the largest 
and most up-to-date catalogue of marine stock assessments available. For each assessed stock 
(comprising a specific species and geographical location), the relevant data consist of yearly 
estimates of the total biomass. These data integrate multiple independent sources of information 
beyond catch data, such as species-specific biological information and the results of research 
surveys, and are consequently regarded as more accurate estimates of population abundance. 
To obtain meaningful statistics and facilitate comparison to results in the main text, we focus 
on those assessments that have at least 30 consecutive years of data and that have total biomass 
estimates measured in units of tonnes. Out of the 331 assessments in the database, 199 satisfy 
these criteria. 

Figure [S6] shows the statistics of the normalised growth rate fluctuations calculated for the 
populations in the RAM database. As shown, the central findings in the main text hold true. 
Namely, the normalised growth rate fluctuations follow a double-exponential (Laplace) distri¬ 
bution with exponent close to y/2, both in the case when data from all assessments are pooled 
together (Fig.|S6)A) and at the individual population level (Fig.|S6ft-C). This analysis provides 


evidence that our results are not artefacts of biased sampling or observational error, but rather 
reflect the true statistical patterns of the underlying population dynamics. 
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Fig. S1: Illustration of non-constant catchability. (A) Relationship between underlying population 
abundance and catch under two different catchability scenarios: hyperstability (blue) and hyperdepletion 
(red), in which catchability increases or decreases at low population abundances, respectively. These 
scenarios should be contrasted with the black curve, representing the assumption in the main text that 
catch is linearly related to abundance. (B) Corresponding catchability functions f(nt) described by 
Eq.|S5|(rcd), Eq.[S7](blue), and the constant catchability relation /(rq ) = a (black). 
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Fig. S2: Validation of the distribution of growth rate fluctuations under non-constant catchability. 

Each curve represents the histogram of aggregated normalized growth fluctuations derived from the 
population abundances {nf 1 }, which are assumed to be related to the catch data {xp } in the LME 
dataset according to the given catchability relation: baseline, hyperstability, or hyperdepletion. 
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Fig. S3: Model validation for non-constant catchability. Counterparts to Fig.[2]4 for the hyperdeple¬ 
tion (left) and hyperstability (right) scenarios. 
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Fig. S4: Statistics of the normalised growth rate fluctuations of individual species. (A and B) 

Examples of the cumulative distributions for (positive fluctuations) and (negative fluctuations, in 
absolute value) along with the corresponding best exponential fits, where the exponent /3 in 1 — e~^t is 
obtained by maximum likelihood estimation. (C) Distributions of the best fit exponents of the individual 
species, for all species in all 66 ecosystems under consideration. The a now indicates the exponent \/2 
obtained for the aggregate of all species (main text, Fig. [Ip), whereas the continuous and dashed lines 
indicate the mean exponent of individual species for positive and negative fluctuations, respectively. 
(D) Distributions of the /;-values calculated using the Kobwogorov-Smirnov test between the empirical 
distributions of and their best exponential fits, represented in panel C. Assuming, as usual, that 
the hypothesis is not rejected for /;-values larger than 0.05, it follows that the normalised growth rate 
fluctuations of most species are consistent with an exponential distribution with exponent close to y/2. 
This confirms that the scaling for the aggregate data remains valid for individual species. We have 
verified that this conclusion also holds using the Akaike Information Criterion (AIC), which is another 
independent test of goodness of fit. According to the AIC, over 70% of species’ growth rate distributions 
more plausibly follow a Laplace distribution compared to the null hypothesis of a normal distribution. 
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Fig. S5: Statistics of non-normalised growth rates, 77. (A) Non-normalised growth rates aggregated 
over all species, where the distributions of positive growth rates (solid symbols) and negative growth 
rates (open symbols) are shown separately. This panel represents the non-normalised counterpart of the 
normalised growth rates considered in Fig. [Ip of the main text. (B) Distribution across all species of each 
species’ standard deviation of the growth rates, < 7 . The distribution of non-normalised growth rates (panel 
A) has a heavier tail than the Laplace distribution representing normalised growth rates; this difference 
is due to heterogeneity across species, which is mainly reflected in the relatively broad distribution of the 
standard deviations (panel B). Thus, the normalisation introduced in our analysis is an important step to 
reveal the scaling behaviour identified in this study. 
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Fig. S6: Statistics of normalised growth rate fluctuations in the RAM Legacy Stock Assessment 
Database. (A to C) Counterparts to Fig. [TJ" (main text), Fig.|S4p, and Fig.|S4p. respectively. 
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Fig. S7: Statistics of population abundances. (A and B) Examples of the cumulative distributions for 
[In xy — [hi xy)] 1 (positive deviation from the average) and [In xy — (hi xy)]' (negative deviation from 
the average, in absolute value) along with the best exponential tits for scaling exponents obtained by 
maximum likelihood estimation. These fits correspond to power-law probability distributions for the 
positive and negative deviations of the abundance itself, P(x) = Pox -1 where Po = j3d i ' ]nx K (C) 
Distributions of the best fit exponents of the individual species, for all species under consideration. In¬ 
set: cumulative distributions of the corresponding /;-values calculated using the Kolmogorov-Smirnov 
test, where the continuous (dashed) line corresponds to positive (negative) deviations. (D) Counter¬ 
part of panel C for best fits of In xy by normal distributions, corresponding to lognormal distributions 

, - , -(In a;-(In a )) 2 

P(x) = (xo\/2tt) e 57 ® for the abundances, where the standard deviations a = cr(may) 
are determined by maximum likelihood estimation. (E) Distributions of the average log-abundance per 
species, (In ay), and of the full time-series of the log-abundances, In ay, aggregated over all species. 
Panel E indicates that the relative species abundances are mainly determined by differences between the 
average populations of different species rather than by their time variations. On the other hand, the distri¬ 
butions of the abundances of individual species over time arc generally no closer to lognormal functions 
(inset of panel D) than they are to power-law functions (inset of panel C). 






































